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Abstract 

A numerical study of domain wall relocation during voltage switching with different ramping 
times is presented for weakly coupled, doped semiconductor superlattices exhibiting multistable 
domain formation in the first plateau of their current-voltage characteristics. Stable self-oscillations 
of the current at the end of stable stationary branches of the current-voltage characteristics have 
been found. These oscillations are due to periodic motion of charge dipoles near the cathode that 
disappear inside the SL, before they can reach the receiving contact. Depending on the dc voltage 
step, the type of multistability between static branches and the duration of voltage switching, 
unusual relocation scenarios are found including changes of the current that follow adiabatically 
the stable I-V branches, different faster episodes involving charge tripoles and dipoles, and even 
small amplitude oscillations of the current near the end of static I-V branches followed by dipole- 
tripole scenarios. 

PACS numbers: 73.63.-b, 05.45.-a, 73.23.-b, 73.40.-c 
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I. INTRODUCTION 



Nonlinear vertical electron transport in semiconductor superlattices (SL) gives rise to 
a rich variety of dynamical phenomena associated with negative differential conductivity 
(NDC) [l|, U S Ql- wide miniband SL, NDC due to Bragg scattering is the origin of 
self-oscillations of the current through a dc voltage biased SL P, 0] • These oscillations are 
due to recycling of charge dipoles as in the Gunn effect of bulk GaAs. In weakly coupled 
SL, NDC due to sequential tunneling between quantum wells may cause either stable self- 
sustained current oscillations mediated by traveling charge monopoles or dipoles 
or a sawtooth multistable current-voltage (I-V) characteristics associated with static field 
domains ^ 11, Q|. The observed behavior depends crucially on the SL configuration 
(widths of wells and barriers, number of SL periods, boundary conditions in the contact 
regions), doping density and voltage bias Q. 

Although there remain important gaps in our theoretical knowledge of nonlinear transport 
in SL l|, IJ, a basic difference between weakly and strongly coupled SL is the type of balance 
equations describing them. Weakly coupled SL are described by spatially discrete balance 
equations whereas spatially continuous balance equations determine the nonlinear behavior 
of strongly coupled SL j^. Traces of spatial discreteness are the sawtooth I-V characteristics 
and the current spikes during self-oscillations and during domain relocation due to voltage 
switching; see the review of theory and experiments in 

Recently, the dynamic response of the current and the field profile to volta ge s witching 
has been investigated both experimentally Ili) 15, 3|l2| and theoretically jl^, 13 • For 
a SL with a multistable I-V characteristics, each branch thereof corresponds to having the 
domain wall separating the low and high field domains of the field profile (which is a charge 
accumulation layer, GAL) placed at a different well of the SL. We are interested in the 
transition from one stable stationary branch of the 1-V characteristics to another due to a 
step in the applied voltage AV = Vf — Vi (Vi and V/ are the initial and final voltage values). 
Bias steps are turned on during a short time interval called ramping time, which can be zero. 
A bias step increasing the applied dc voltage is referred to as an up jump {AV > 0), while 
a bias step decreasing the applied dc voltage is called a down jump {AV < 0). Bias steps 
contrast with voltage up-sweeps and down-sweeps for which the ramping time is infinitely 
long and the bias increase or decrease is adiabatic. For down jumps, the relocation process 
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of the domain boundary proceeds via a direct motion of the CAL in the direction of electron 
flow. This behavior is confirmed by single-shot time traces of the current response: for values 
of Vf away from regions of bistability, there is an initial displacement current spike, after 
which the current rapidly switches to the stable value. Furthermore, when Vf is near to the 
bistable region, there is an additional intermediate period, in which the current fluctuates 
about a metastable value for a stochastically varying delay time r^, before rapidly switching 



to the stable value in a time r. 



15l |. However, for up jumps, the charge monopole at the 



domain boundary would have to move against the electron flow, and this is only possible for 
small-amplitude up jumps. For larger up jumps, the more complex dipole-tripole scenario 
occurs: the CAL moves one well against the electron flow, then a charge dipole comprising 
one CAL and one charge depletion layer (CDL) is formed at the cathode and, together with 
the old CAL, it moves with the electron flow. The resulting charge tripole exists until the 
old CAL reaches the anode and disappears, leaving only the charge dipole. The CDL of 
this dipole reaches the anode while its CAL moves until its final position at the destination 
branch of the I-V characteristics [la, Il8|. For Vf near the end of a branch, there are 

nil 

pronounced stochastic effects due to shot noise |15L I17|. 

In this paper, we carry out an extensive numerical study of the dynamical response to 
voltage switching and unveil unexpected behavior. We use a discrete sequential tunneling 

n 

model whose detailed description can be found in Appendix A of [2] . Stochastic effects due 
to shot noise will be ignored. We find stable self-oscillations of the current at the end of 
stable stationary branches of the I-V characteristics. These oscillations are due to periodic 
motion of charge dipoles near the cathode that disappear inside the SL, before they can reach 
the receiving contact. We also want to understand how the dynamical response to voltage 
switching is affected by the number of multistable branches of the I-V characteristics, their 
extension and the ramping time necessary to change voltage (from Vi to Vf). Among our 
results, we find a different tripole- dipole scenario than that reported by Amann et al Q| 
(the first phase of the scenario is different). We also find that the ramping time selects the 
tripole-dipole scenario for large up jumps. Suppose that there are several branches of the 
I-V characteristics in the interval between Vi and Vf (large voltage switching). Then there 
are two critical ramping times r^i and Tc2, Tc2 < Td, whose precise values depend on the 
SL parameters in Table H] and on Vi and Vf. For the parameters used in our simulations, 
the critical ramping times are between 10 and 30 /is. If the ramping time is larger than 
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Tel, the current follows adiabatically the stable branches until their end, falls to the next 
stable branch and repeats this process until Vf is reached. For ramping times between Tc2 
and Tci, adiabatic motion over a stable branch is followed by a tripole-dipole scenario until 
the following stable branch is reached. Depending on the number of multistable branches, 
sometimes a stable branch is skipped in this process. Lastly, if the ramping time is shorter 
than Tc2, the final stable branch is reached after only one tripole-dipole scenario even for 
large voltage steps. 

The rest of the paper is as follows. The model we use and details of its numerical integra- 
tion are described in Sections |n] and IIIH respectively. Section IIVI contains the multistable 
I-V characteristics of a SL with realistic configuration and doping density parameters ^. 
We show that the width of the multistability regions increases with voltage while the slope 
of the branches (which is the positive differential conductivity or, in short, PDC) decreases. 
These features affect substantially the dynamic response to switching described in Sections 
fVl and IVIl Finally, the main results of the paper are summarized in Section I VI II 



II. SEQUENTIAL TUNNELING DISCRETE MODEL 

In weakly coupled SL, typically the scattering times (about 0.1 ps) are much shorter 
than the escape times from quantum wells (about 0.01 ns for a SL with a 0.1 meV miniband 
width). In their turn, the latter are shorter than typical dielectric relaxation times (on the 
nanosecond time scale) This implies that the dominant mechanism of vertical charge 
transport is sequential resonant tunneling and that the tunneling current across barriers 
can be considered to be stationary on the time scale of dielectric relaxation. Nonlinear 
stationary and oscillatory phenomena occurring for voltages in the first plateau of weakly 
coupled doped SL have been well described by the spatially discrete model equations (with 
backward finite differences) introduced in [2^ with constitutive relations between sequential 



tunneling current, electron densities and electric field of the type calculated in j2l| using 
stationary nonequilibrium Green functions or in references cited therein) approxi- 

mating Transfer Hamiltonian formulas. See the review [4] for a recent description and further 
justification. The model equations consist of the Poisson and charge continuity equations 
for the two-dimensional (2D) electron density Ui and average electric field —Fi at the ith SL 
period (which starts at the right end of the {i — l)th barrier and finishes at the right end of 
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the ith barrier) 
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Here No, £, — e and eJj^j+i are the 2D doping density at the zth well, the average permittiv- 
ity, the electron charge and the tunnelling current density across the ith barrier, respectively. 
The width of a SL period is / = d + w, where d and w are the barrier and well widths, respec- 
tively. Time-differencing Eq. and inserting the result in Eq. (j2)), we obtain the following 
form of Ampere's law, 

dF 

e^ + J^^i+i = J{t), (3) 
which may be solved with the bias condition for the applied voltage V{t): 
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V{t) 



(4) 



{N + l)f 

The space-independent unknown function J(t) is the total current density through the SL. 
The 2N + 2 independent equations of the discrete model are (0) for i = 1, . . . ,N, (jHl) for 
i = 0, . . . , N and (jH) for the 2N + 2 unknowns nj, Fi {i = 1, . . . , N), Fq and J, provided we 
have + 1 constitutive relations linking the tunneling current (z = 0, . . . , A^) to the 

electron densities and electric fields. To calculate the tunneling currents across SL barriers, 
we use explicit formulas provided by the Transfer Hamiltonian method when the scattering 
broadening is much smaller than the subband energies and chemical potentials j^: 
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As boundary tunne 
and field as in Ref. 



ing currents for i = and A^, we adopt linear relations between current 



18l |. In these formulas, a is the contact conductivity (assumed to be the 



same at both contacts for simplicity), m* the effective mass, T the temperature, fc^ and h 
are the Boltzmann and Planck constants respectively, and the "forward tunneling velocity" 
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v^f^ is a sum of Lorentzians centered at the resonant field values Fcu = {^cu — £ci)/{.d)\ 
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(10) 
(11) 
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In these equations, Cj indicates the jth subband in a well, Ec^ is its energy measured from the 
bottom of the well, is the scattering width, 7^ is a dimensionless transmission probability 
across the ith barrier, and eVf, is the barrier height in absence of potential drops. Typical 
values of these parameters are shown in Table HI 



TABLE I: Parameters of the 9/4 SL in Ref 



N Nd dw/ds 7 ^Ci ^Cs Vb 

(cm-2) (nm/nm) (meV) (lO^^a Kg) (meV) (meV) (meV) (V) 



40 1.5 X 10" 9.0/4.0 8 



8.43 



44 



180 410 0.982 



To carry out numerical integrations of the discrete model, the explicit formulas ^ - ()14|) 
are much better than numerically calculated tunneling currents such as those obtained in 



21| from Green function calculations. Furthermore, explicit tunneling currents are better 



suited for analysis of the discrete model equations. These reasons to favor the previous 
explicit formulas for the sequential tunneling current are not offset by claims that one type 
of derivation (Transfer Hamiltonian or Green functions) agrees better with first principles: 
both derivations involve similarly drastic simplifications and the resulting formulas agree 
similarly well with experiments. It turns out that the type of solutions of the discrete model 
depends on the qualitative features of Jj^j+i as a function of Fj, rii and nj+i, not on detailed 



quantitative features. The formulas for Ji^i+i obtained using different derivation methods 
yield a tunneling current similar to that in Fig. ^ which is why we obtain similar results. 
See discussions in Ref. m. 



III. NONDIMENSIONALIZATION AND NUMERICAL INTEGRATION 



For numerical treatment, it is convenient to render the equations dimensionless. We have 
used the following definitions 

t Jut 
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The values Fm and Jm are defined as the field and current density at which the tunneling 
current Jj_»i+i of (0) reaches its first relative maximum, provided rii = nj+i = A^'^. With 
these definitions, the model equations are 
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(22) 



sFm ttH^Nd ksT' 

are dimensionless parameters. Their values for the SL described in Table H] are given in 
Table IHI Figure ^ shows the tunneling current density as a function of a homogeneous field 
profile Fi = F when the electron densities are set equal to the doping density in all the SL 
wells. We observe several relative maxima, the first of which yields the values of Fm and 
Jm- The field intervals between two consecutive maxima roughly correspond to plateaus in 
the I-V characteristics. 
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TABLE II: Typical scales for T = 5 K. 



Fm Jm vm xo to V po Pc Vb 

(kV/cm) (A/cm^) (m/s) (nm) (ns) (-) (-) {VL m) (V) 



Jul eFul eFm eNo m*kBT IFm ^ 



eNo eNo Jm sFm ttTi Nd svmNd 



3.945 3.127 1.691 2.494 2.066 5.212 0.111 12.62 0.205 





100 200 300 

(kV/cm) 



400 



FIG. 1: Tunneling current density for rij = rij+i = as a function of field Fi = F, showing that 
Fm « 3.945 kV/cm and Jm » 3.1269 A/cm^ for T = 5 K and the SL values of Table U 

In order to numerically solve the discrete model, we first sum |TBj) from z = to and 
use the bias condition (fTHj) to calculate J'. The result is 



J 



d(t){t) 
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dt N + 
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To solve (|16| ) - (j2H) together with the initial condition 
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is equivalent to solving (fT7|) - (j^H) plus the following equation instead of (fT^ : 
dTi dcj) I ^ ^ N ^ .... 

The new system of equations also satisfies the bias condition psp . This can be checked 
by adding all equations ()25|) from i = to which implies that Y.iLo-^i ~ + 
is a constant, equal to zero because of the initial conditions. To solve our dimensionless 
equations, we have used an embedded Runge-Kutta method of order 7(8) with step-size 
control and error estimate, checking the results independently by means of an implicit BDF 
(backward differentiation formula) method of order 1 to 4, solved by means of Newton- 
lap hson iterations. These methods are more accurate than those used by Amann et al 



IV. I-V CHARACTERISTICS 



We have constructed numerically the first plateau of the I-V characteristics of the SL 
whose parameters are compiled in Table |l] at a temperature of 5K and a contact resistivity 
1/a = 25.2 Qm {a = 0.5). We depict in Fig. El the portion of the stable static branches 
obtained by voltage up-sweep from zero volts to the end of the first plateau and also the 
portion obtained by voltage down-sweep from high voltage values. Note that the upper 
parts of odd numbered branches (from the 23rd to the 33rd) do not appear because they 
correspond to regions of bistability and are skipped during voltage up-sweep. Likewise, the 
central parts of branches 35th to 39th are not shown because they correspond to regions 
of tristability and are skipped during voltage up-sweep and down-sweep. All these curves 
can be entirely shown by up-sweeping starting from the parts thereof shown in the figure. 
The voltage distribution for different locations of the domain boundary for a fixed applied 
voltage in regions of multistability is shown in Figures 8 and 9 of the review paper |^, which 
contains further background discussions and references to pertinent literature. The basin 
of attraction of the static branches shown in Fig. |21 is very small near their ends and it is 
mainly in these regions that our high precision numerical methods make the difference with 
previous calculations; see Fig. El We have found that the I-V branches may undergo Hopf 
bifurcations to small-amplitude oscillatory solutions near their upper ends, as we explain in 
Section IVIl 
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FIG. 2: (Color online) I-V characteristics of the 40-period 9-4 SL of Ref. 8] obtained by up- and 
down-sweeping adiabatic processes for V G [0,4] V. Parameters correspond to Table at 5K and 
with cathode resistivity of 25.2 Jim (a = 0.5). The branch number increases with voltage: the ith 
branch has a CAL separating low and high field domains which is located at the {N — i -|- l)th well. 



According to their stability feature, we can distinguish four different types of I-V branch: 

1. The first branch Bi is singly stable from V = until its upper part coincides with the 
lower end of branch i?2, giving rise to a narrow region of bistability. 

2. From branch B2 to i?2i, we observe that the branch length increases (the extension of 
B2 is 0.122 V whereas that of -B19 is 0.177 V) while their slope (the PDC) decreases. 
We see that the low voltage branches have a steep slope similar to that below the 
first maximum of the homogeneous tunneling current density in Fig. ^ whereas the 
high voltage branches have a smaller slope similar to that of the low part of the 
ascent to the second maximum in Fig. ^ Why? The field profile in the low voltage 
branches is similar to that of a solution with spatially uniform field, which therefore 
obeys Ji^i+i{F, No, Nd) = J, thereby corresponding to the first branch of the curve in 
Fig. n] Similarly, high voltage branches are close to spatially homogeneous field profiles 
satisfying the same relation, but now the field profile corresponds to the third branch 
of Fig. ^ Branches with intermediate voltages are a combination of low and high field 
domains, and therefore their slopes are interpolations between low and high PDC. The 
central part of each branch is singly stable, while their two ends are bistable. As the 
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FIG. 3: (Color online) Bistability between Branches 5 and 6 of the I-V characteristics in Fig. |2l 
{V £ [0.55,0.6] V). The stationary solution becomes unstable to a small-amplitude oscillatory 
solution at the upper end of Branch 5. 



branch number increases, the central part shrinks and the bistable regions including 
the ends of the branch grow. 

3. Branch B22 is the last one having a central singly stable region and the first one with 
a tristability region near its upper end. In the tristability region, branch B22 coexists 
with the central part of branch B23 and the lower part of branch B24,. 

4. Branch B23 is bistable except in its central part and its upper end where it is tristable. 
The lower part of this branch is still bistable. 

5. Branches B2A to B^g are tristable in their central parts and ends, but they have two 
bistability regions. As voltage increases, the tristability regions grow at the expense 
of the bistability regions. 

6. The last branch of the first plateau, B^q, has a tristable lower end, a bistable central 
region and it is singly stable from there until it reaches the second plateau in its upper 
end. 
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V. LARGE SWITCHING: EVOLUTION OF J{V) ALONG THE I-V CURVE 



In this Section, we describe the dynamical response of the SL to a voltage switching 
V{t) = Vi + AVtH^Tr - t)/Tr + AV H (t - Tr), m which AV is constant, and H{t) = 1 if 
t > 0, H{t) = if t < is the Heaviside unit step function. We select the initial voltage 
in the central part of one branch and the final voltage Vf = Vi + AV in different parts of 
another stable branch, so that several branches can be found between Vi and Vf. We have 
found different scenarios. 

A. Switching from bistable branches: modified tripole-dipole scenario 

Let us choose Vi = 0.83V, in the central singly stable part of branch Bg and Vf = 1.37V, 
on the central part of branch i?i4. These branches are bistable but they have a central part 
for whose voltages no other static solution is stable (the branch is singly stable there). Fig.|3] 
shows the dynamical response of the total current density J{t) to voltage switching with two 
different ramping times. For large enough ramping times (not shown), the current density 
J{t) follows adiabatically the I-V characteristics. Below a first critical value of the ramping 
time, J{t) cannot reach the upper end of the branches and it falls to the lower part of 
the following branch in the voltage range where both branches are stable. For each branch 
crossing, this fall occurs via the following modified tripole-dipole scenario whose current 
trace is depicted in Fig. 13 

Phase 1 The CAL separating low and high field domains stay in the same well while the current 
density increases until it surpasses the critical value at which J = crF intersects the 
curve in Fig.[T[ J = Ji^i^i{F, Nd, Njj) on the second branch thereof (Jc ~ 2.9 A/cm^). 

Phase 2 A charge dipole wave is created at the cathode and it starts moving to the anode while 
the old CAL also moves towards the anode. After a short transient during which the 
current decreases, the whole structure, a charge tripole moves rigidly towards the 
anode. 

Phase 3 The old CAL reaches the anode leaving a charge dipole moving (after a short transient 
during which the current further decreases) rigidly towards the anode at a lower speed 
than the tripole. 
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FIG. 4: (Color online) Total current density during voltage switching from Vi = 0.83V to Vf = 
1.37V (seven branches) for ramping times: (a) r,. = 30;us (dimensionless value: ?> = 2 x 10^), and 
(b) Tj. = 15/Us {fr = 10^). Thick black line: upper part of the I-V branches, red thin line: lower 
part of the static I-V branches, thin blue line: response curve {V{t),J{t)). In (a), the current 
follows adiabatically the I-V curve, whereas in (b), the ramping time is so short that the 4*'' and 
6*^ branches are skipped. The dimensionless cathode conductivity is fi = 0.5. 

Phase 4 The front part of the dipole, a CDL, reaches the anode leaving only the CAL in its 
back part. The current density increases to higher values corresponding to the next 
stable static solution while the CAL moves towards its final destination separating low 
and high field domains of the stable static field profile. 

Phase 2 is characterized by double-peaked current spikes (corresponding to the well-to-well 
jumps of the two CAL) while Phases 3 and 4 exhibit single current spikes; see Fig. El Phases 
2, 3 and 4 are exactly as described in [18] (who did not describe Phase 4) and corrected in 
[16[ (who added Phase 4). In these previous works, Phase 1 was characterized by a one- well 
motion of the CAL towards the cathode. Note that experimental observations refer only 
to the behavior of the current density, not to the motion of domain walls, and therefore 
they cannot discriminate between our Phase 1 and the different one reported in To 
explain the differences in Phase 1, we note that, according to Fig. 3 of Ref. [2^, a CAL 
may move against the electron flow if the current is large enough and the doping density 
is sufficiently large. However, the critical current for this motion would be exponentially 
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FIG. 5: Time trace of the total current density J{t) after sudden voltage switching characterizing 
the tripole-dipole scenario. Here Vi = 1.0 V (for t < 0) and AV = 0.2 V. 

close to the maximum current for which CAL exist if the doping density is even larger, 
thereby eliminating in practice the possibility for CAL to move against the electron flow 
under current bias. In the simulations by Amann et al, the interval of currents allowing 
CAL motion against the electron flow was relatively wide (cf. their Figure 4), whereas we 
found that it was negligible for the tunneling current and doping density used in the present 
calculations. Then the increase of voltage due to switching is compensated by simply a 
current increase in Phase 1, without CAL motion, as indicated in Fig. El When the current 
surpasses its critical value (see the current spike over the maximum current value for static 
branches in Fig.|^, Phase 2 begins: a dipole is shed from the cathode and it moves on towards 
the anode tog ether with the old CAL. This and subsequent Phases are as in the previous 
works The behavior of the current can be explained using singular perturbation 

ideas as described in the reviews Since the final position of the CAL approaches the 

cathode as the voltage increases, the dipole trip becomes shorter as shown in Fig. |^ 

We have seen that switching from to B14 occurs by a succession of tripole-dipole 
scenarios. What happens if we continue decreasing the ramping time? It turns out that we 
start skipping branches. Note that the tripole-dipole process is quite fast in Fig. |H^a) and 
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FIG. 6: Field profiles corresponding to Fig. ^during relocation of the domain wall, (a) Detail of 
the emission and the travel of the CAL during the second relocation stage in Fig. El^a). (b) The 4*^* 
branch in Fig. IHb) has been skipped, therefore the 4*^ CAL is missing and there is no relocation 
stage. 

it ends at a voltage smaller than the bistability interval near the end of the corresponding 
I-V branch. At that voltage, only one branch is stable and J{t) follows this static branch 
until the next dipole emission. However, recall that the bistability intervals grow at the 
expense of the singly stable central part of the branches as their voltage increases. When 
the ramping time is decreased below a critical value (which depends on the branch), the 
voltage at the end of a process of dipole emission and travel may be in the bistability range 
of two branches. The CAL of the dipole then stops at the position corresponding to the 
static branch with lower current which is closer to the cathode than the CAL of the branch 
with higher current. This is observed in Fig. Efb). Note that the branches Bu and i?i3 have 
been skipped. Fig. Efb) shows that the corresponding tripole-dipole process is longer when 
one branch has been skipped. 

B. Same switching range for bistable and tristable branches 

It should be clear from the previous discussion that the intervals of bistability have 
great influence on the dynamic response to voltage switching. To make this clearer, we 
have depicted in Figures Et^a) to (h) the dynamic response to a voltage switching of width 
AV = 0.5 V with ramping time Tj. = 30/is (as in Fig. EJ^a)) for different initial voltages 
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Vi. We observe that the tripole-dipole scenario occurs for small Vi, branches start to be 
skipped as Vi increases, and the tripole-dipole process disappears at even larger voltages at 
which the tristability range of the branches is very large. Of course, the occurrence of the 
tripole-dipole scenario depends on the ramping time: it reappears again as the ramping time 
decreases sufficiently. 
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FIG. 7: (Color online) Voltage switching with = 0.5 V and ramping time = 30 for 
different Vi on the I-V characteristics: (a) Vi = 0.5 V, (b) = 1 V, (c) Vi = 1.5 V, (d) Vi = 2V, 
(e) Vi = 2.5 V, (f) = 3 V, (g) V = 3.5 V, (h) = 4 V, AF = 1.5 V. In (a), (b), (g) and (h), 
the current follows the I-V curve along the upper part of each branch, whereas J ranges from 6 to 
16 A/cm^. In (c), a branch (the 5*^) is skipped for the first time. In (d) and (e), branches 3 and 
5 are skipped. In (f), the upper part of the branch is reached but the G*'^ is skipped. 
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VI. VOLTAGE SWITCHING TO Vf NEAR THE END OF A BRANCH 



In this Section, we report the current response to voltage switching from Vi to a voltage 
Vf close to the end of the same branch. Firstly, we shall select branches near the end of the 
plateau, which were the only ones considered in Ref . 3] (see Section VI) . Secondly, we shall 
select branches near the beginning of the plateau and observe a rather different behavior. 



A. High voltages near the end of the I V plateau 

Fig. ISl depicts the current response to a voltage switching starting with a in the bistable 
part of branch S37, closer to its end. Kee ping a ramping time of 100 ns, we observe similar 
behavior to that reported by Amann et al lo\: (i) if Vf < Vth (the end of the static branch), 
the current remains on the same branch but the time it takes to settle in its final value 
increases as Vf approaches Vth'-, (ii) if Vf > Vth, the final state is on the next branch, and 
the transient stage lasts longer as we approach Vth, cf. Fig. IHl^a). Similarly, the longer the 
ramping time is, the longer the transient before the current drops to that of the following I-V 
branch seems to be, as indicated in Fig.E^b). This Figure shows the influence of the ramping 
time on the current response to voltage switching with Vf close to Vth for I-V branches near 
the end of the plateau. We have selected now Branch ^35 and changed the ramping time 
from small to large for two different Vf close to Vth, one larger than Vth, the other smaller. 
For Vf > Vth, "we observe that the current eventually drops to the lower value on branch 
i?36, but the transient stage lasts longer as the ramping time increases. For V/ ~ K < Vth, 
the basin of attraction of B^^ is so small that the current eventually drops to its final value 
on Branch B^q. However, the way in which this happens depends strongly on the ramping 
time: (i) if the ramping time is too small, the current drops rapidly and the final state is on 
Branch B^q is reached soon; (ii) for intermediate ramping times, the current oscillates about 
the static Branch B35 before it falls to Branch B^^, and (iii) for large ramping times, the 
current seems to settle down to the static value on Branch B35 and it resists much longer 
before it eventually drops down to Branch B^q. 

A distinct feature of the current response to voltage switching with Vf near the end of a 
static branch is that, in case (i), the final stable state for Vf < Vth very close to Vth, may 
be oscillatory, not stationary. This is suggested by the oscillations in Figs. |H1 and El and 
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FIG. 8: Current response to voltage switching from Vi = 3.68 V to five different values of Vf near 
Fth (1: 3.74, 2: 3.747, 3: 3.7472732, 4: 3.7474 and 5: 3.76 V). The ramping time is = 100 ns and 
we have set t = at the end of voltage switching. Note the oscillatory behavior for Vf = 3.747 V. 



further confirmed by the linear stability analysis of Appendix 1X1 There it is shown that the 
static branch loses stability at some Vo < Vth because the real part of two complex conjugate 
eigenvalues becomes positive for Vo < V < Vth- This is clearly seen in Fig. [TUl which depicts 
the eigenvalues of the linear stability problem about the static branch B^^ for 20 different 
voltage values close to Vth- For this branch, Vo ~ 3.5547 V. 



B. Current response to switching near the end of low voltage branches 

Fig. shows that switching near the end of low voltage I-V branches is more complex 
than that described previously. It turns out that the current drop to the next branch 
may occur via the tripole-dipole scenario, unlike in the numerical simulations by Amann 
et al (who always selected high voltage static branches having intervals of tristability), 
but according to the experimental observations by Rogozia et al (cf. Fig. 9 of Ref. 16[)- 
Another discrepancy between the numerical simulations of Ref. and experiments is that 
the current spike accompanying dipole emission is much taller in the simulations (twice 
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FIG. 9: Current response to voltage switching for Vi = 3.5 V, two different values of Vf and three 
different ramping times in each case: (a) Vf = 3.555 V > Vth, for Tr = 1 ns, 1 ^s and 100 ^s. (b) 
Vf = 3.5545 V < Vth; for = 0.1 ^s, 1 ^s and 100 ^s. In all cases, we set t = at the end of 
voltage switching. For = 100 /is, the current undergoes a small-amplitude oscillation in the case 
of Fig. Mh). 
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the maximum current of the static branches instead of the experimentally observed 20% 
increase). In our simulations, the current spike accompanying dipole emission is much 
smaller than in the previous calculations by Amann et al, but this is due to the different 
cathode conductivity and critical current for dipole creation. 

The appearance of the tripole-dipole scenario during voltage switching occurs up to volt- 
ages corresponding to branches starting to display tristability. If we fix Vf sufficiently close 
to Vth, Vf < Vfh, and change the ramping time, we observe a peculiar behavior. For all 
ramping times, the relocation of the domain wall separating the low and high field parts of 
the field profile happens via the tripole-dipole scenario, as shown in Fig. IT^ a). However the 
time trf at which the tripole-dipole scenario starts is not a monotone function of the ramping 
time: it seems that td may have local maxima and minima as a function of r^. We have 
observed that increases with r,. up to ~ 2.2/is. Then decreases with up to at least 
7 yus. Then td increases again for larg Gr Tj-, clS indicated in Fig. I12r a). We have checked that 
td is again smaller for a ramping time of 100 /iS than for r^. = 30/iS. 
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FIG. 10: Evolution of the 41 eigenvalues determining linear stability of the static branch S35 for 
20 values of the applied voltage from V = 3.54 to 3.56 V. 

To ascertain the origin of the non-monotonic behavior of as a function of the ramping 
time, we observe that, as the voltage increases with time, the current becomes oscillatory 
before the tripole-dipole scenario begins. The shape of the current oscillations and their 
local period also change as the time elapses, which is clearly seen in Figs. II II and fTW b) . The 
latter figure also shows that, for similar ramping times, the current may drop to the lower 
I-V branch or continue oscillating for a longer time. Fig. 1131 depicts the field profiles during 
the oscillations of the current shown in Fig. IT^ b) for Vf < Vth- We see that the current 
oscillations correspond to the periodic formation of a small field pulse at the cathode and its 
advance towards the anode over a few SL periods before it shrinks and vanishes. Eventually 
as the voltage increases with time, a dipole succeeds in growing sufficiently to detach itself 
from the cathode region and trigger a tripole-dipole event, bringing down the current to its 
stable value in the next I-V branch. Fig. shows that the same mechanism is responsible 
for similar small amplitude current oscillations for Vf > Vth- 

Another clue to the different current response to switching at low voltage values is offered 
by linear stability analysis of the I-V static solution branches. At low voltages, they have 
a stationary instability at a voltage smaller than Vth- At that voltage, a stationary branch 
bifurcates and this secondary branch undergoes a Hopf bifurcation to a small amplitude 
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FIG. 11: (Color online) Current response to voltage switching between Branches 5 and 6 of the 
I-V characteristics in Fig. IJJfor Vi = 0.56 V, ramping time Tr = 10 ^s, and four different values of 
Vf near Vth: 1: 0.59, 2: 0.591, 3: 0.59203 and 4: 0.593 V. Note that curve 3 is shifted down 0.2 
A/cm^ for the sake of clarity. 

oscillation. The field profiles for these oscillatory solution correspond to the periodic emission 
and motion of a small field pulse (charge dipole) confined to the region near the cathode, 
which is very similar to the confined Gunn effect in bulk n-GaAs and ultrapure p-Ge j^]. 
Figures ITM a) and (b) show that the eigenvalue determining the linear stability of the static 
branch is real and it vanishes at a certain voltage smaller than Vth- The other eigenvalues 
with large real parts are complex and have negative real parts. The gap in Fig. [T^T b) 
corresponds to the oscillatory instability observed in the numerical simulations of the discrete 
model. It is plausible that the stationary branch that bifurcated from has a secondary 
oscillatory instability at that voltage, but a more detailed study is necessary before this 
can be ascertained. Thus we may have the following succession of bifurcations from the 
static branch Br, as the voltage increases towards Vth'- (i) stable -B5, (ii) small amplitude 
static branch issuing from B^, (iii) Hopf bifurcation from the bifurcating static branch, (iv) 
annihilation of the oscillation before or at Vth (the end of -B5). During voltage sweeping. 
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FIG. 12: (Color online) Current response to voltage switching between Branches 5 and 6 of the 
I-V characteristics in Fig. [21 for Vi = 0.56 V. (a) Vf = 0.593 V > Vth and seven different ramping 
times from 1 to 25 /is. (b) Details of current response for Vj = 0.59203 V < Vth and three short 
ramping times: 0.1, 2.1994 and 2.1995 fis. 



the current should go through this succession of bifurcations and there is no reason why the 
time td should be a monotone function of the ramping time Tr- 



VII. CONCLUSION 



The relocation of the domain boundary in weakly coupled doped SL is substantially 
affected by the ramping time over which the voltage is switched and by multistability of the 
initial and final static I-V branches involved in switching. Let us consider voltage switching 
leaving several I-V branches between the initial and final voltages. If the ramping time is 
very long, the current simply follows adiabatically the change in voltage during switching, 
much as in up and down voltage sweeping. If the ramping time is very short, each branch 
jump during switching is achieved by a modified tripole-dipole scenario: a CDL is formed at 
the cathode, it moves towards the anode producing a second CAL behind it. Together with 
the old CAL, the resulting charge tripole moves towards the anode until the first CAL and 
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FIG. 13: Evolution of the electric field profile corresponding to Fig. I12r b'). where Vi = 0.56 V 
and Vf = 0.59203 V < Vth, for the two (very similar) ramping times (a) r,. = 2.1994 ^us and (b) 
Tr = 2.1995 ;us. In (a), the CAL is emitted after a short oscillatory transient. In (b), a CAL is 
emitted from the cathode, it disappears in the interior of the sample and this behavior is repeated 
periodically. This is similar to Gunn effect oscillations confined to one part of the SL. 



the CDL reach it. Then the remaining CAL moves to its final position corresponding to the 
new static I-V branch. For intermediate ramping time, and provided the I-V branches have 
wide intervals of bistability, the tripole-dipole scenario may be skipped, thereby occurring 
every other branch jump. If the final voltage after switching is very close to the end of a 
I-V branch, the current eventually drops to its value at the following static branch, but it 
can remain a long time on the initial static branch (or it oscillates about it in case there 
is an oscillatory instability) if the ramping time is sufficiently long. The time at which 
the tripole-dipole scenario begins and the current drops to its stable value at the next I-V 
branch is not a monotone function of the ramping time. 
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FIG. 14: Evolution of the electric field profile corresponding to curve number 6 in Fig.^J^a), which 
has longest oscillatory interval before a C AL is emitted from the cathode (ramping time r,. = 20/is) . 
(a) Detail of the oscillatory transient regime, (b) Detail of the profile near the cathode (closest 15 
wells) when the CAL is finally emitted. 
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FIG. 15: (a) Evolution of the 41 eigenvalues determining linear stability of the static Branch 
for 20 values of the applied voltage from V = 0.56 to 0.6 V. The insets show the overall picture 
with all the eigenvalues and a zoom of the region near the imaginary axis, (b) Real part of the 
eigenvalues as a function of voltage. The gap in the figure corresponds to the oscillatory instability. 

APPENDIX A: LINEAR STABILITY OF THE STATIC I V BRANCHES 

Let ({jF*}^g, J*) a stationary solution of (fTI)jl - (PT|) under dc voltage bias (f){i) = 0, Vt. 
In these equations, we shall eliminate the electron density in favor of the field by using the 
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Poisson equation. Then the tunnehng current is a function of the electric field profile such 
that 



J^^^+l{{:Fn) = J\ 2 = 0,...,iV, 



N 



^j;* = (iV + i)0. 



(Al) 
(A2) 



i=0 



Let ({/i(t)}i=0'i(^)) ^ disturbance from the static solution: 



Ti{i)=T: + efi{i), J{i)=J* + ej{i). 



(A3) 



Then, the linear equations about the static field profile and the static current density are 



dJi 



fi + l 



0,...,Ar, (A4) 



(*) 



up to 0(e) terms. We have set (*) = {T* , . 

Let us now assume fiit) = e^^fi, j(t) = e^^ j. Then we obtain 



A/i =j - ft 



n+l 



9^ 



-f. 



(*) 



(A5) 



which can be written in matrix form as Af = j — A ■ f , with ai^i = 



(*)■ 



A 



Jo J 

h 



\fN ) 



J 



(^0,0 (^0,1 
^1,0 0,1,1 0,12 
0-2,1 02,2 



V 



on-i,n 

On,N-1 On^n J 



( f \ 

Jo 

fl 

\fN) 



(A6) 



The boundary conditions for z = and i = N yield ao,o = cr, Oq i = 0, aN,N-i = 
and otv.at = o"(n|^ + J-'^/u). On the other hand, the bias condition (fT8|) becomes 



-a 



N 



(A7) 



i=0 



Then Af + A ■ f = j implies f = (AI + A)"^ • j. (jAl| means 
the vector f is zero, therefore we have 



5:(AI + A)-i-l = 0, 



that the sum of the entries of 



(A8) 
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because j ^ 0. The left hand side of ()A8|) is polynomial of degree in A, having therefore 
N zeros. For computational purposes, it is better to rewrite the system flA6|) - ()A7|) in the 
form Af = B • f . This can be achieved adding the rows in ()A6|) and using ()A7|) : 



N 



>^J2fi = i^ + ^)j~ ("0,0 + ai,o)/o - (ao,i + + ct2,i)/i - • • ■ 

... — (oj-i^i + Qi^i + ai+i^i)fi — ... — {aN-i,N + aN,N)fN = 0. 



i=0 



(A9) 
(AlO) 



Defining Si = ai_i,i + ai,i + Oi+ii, i = 1, . . . , N - 1, sq = ao,o + ai,o and sn = aN-i,N + aN,N, 
we have 

1 



N 



and therefore 





J = 


(A 




j 


1 

N+1 


[^1 





N + 



(All) 



Sq Si S2 

So Si S2 



Sn 
Sn 



V 



/ f \ 

Jo 

fl 



[fNj 



(A12) 



So Si S2 . . . Sn 

which is of the form j = jy+T^ ' ^- Substituting this expression in ()A6|) . we obtain a matrix 
equation of the type Af = B ■ f , namely 



A 



Jo 

fl 

[In) 



1 



N+1 



So Si S2 . . . Sn 

So Si S2 . . . Sn 

\So Si S2 ■ ■ ■ Sn J 



\ 



0-0,0 0-0,1 

0-1,0 Oi^i ai^2 








on-i,n 
on,n~i on,n j 



( f \ 

JO 

{In) 



The matrix B is equal to the matrix S/(A^ + 1), except in its three main diagonals, where 
bi,j = + 1) - Oij: 



A 



fo 
fl 



bo,0 ^0,1 ^0,2 • • • bo,N 
bl,0 ^1,1 &1,2 • • • bi^N 



\ f f \ 

Jo 

fl 



(AM) 



^fNJ \ bN,0 bN,l bN,2 ■ ■ ■ bN,N J \ f^ j 

The (A^ + 1) X (A^+1) matrix B has a zero eigenvalue (add its rows), and its other eigenvalues 
are the zeros of the polynomial ()A8|) . They have been depicted in Fig. ITUl 



[*] E-address bonillaOing. uc3m.es. 



26 



[1] A. Wacker, Phys. Rep. 357, 1 (2002). 

[2] L. L. Bonilla, J. Phys.: Cond. Matter 14, 341 (2002) 

[3] G. Platero and R. Aguado, Phys. Rep. 395, 1 (2004). 

[4] L. L. Bonilla and H. T. Grahn, Rep. Prog. Phys. 68, 577 (2005). 

[5] M. Biittiker and H. Thomas, Phys. Rev. Lett. 38, 78 (1977). 

[6] K. Hofbeck, J. Grenzer, E. Schomburg, A. A. Ignatov, K. F. Renk, D. G. Pavel'ev, Y. Koschuri- 

nov, B. Melzer, S. Ivanov, S. Schaposchnikov, and P. S. Kop'ev, Phys. Lett. A 218, 349 (1996). 
[7] J. Kastrup, R. Klann, H. T. Grahn, K. Ploog, L. L. Bonilla, J. Galan, M. Kindelan, M. 

Moscoso, and R. Merlin, Phys. Rev. B 52, 13761 (1995). 
[8] J. Kastrup, R. Hey, K. H. Ploog, H. T. Grahn, L. L. Bonilla, M. Kindelan, M. Moscoso, A. 

Wacker, and J. Galan, Phys. Rev. B 55, 2476 (1997). 
[9] D. Sanchez, M. Moscoso, L. L. Bonilla, G. Platero and R. Aguado, Phys. Rev. B 60, 4489 

(1999). 

[10] L. Esaki and L. L. Chang, Phys. Rev. Lett. 33, 495 (1974). 

[11] Y. Kawamura, K. Wakita, H. Asahi, and K. Kurumada, Jpn. J. Appl. Phys. 25, L928 (1986). 
[12] H. T. Grahn, K. von Klitzing, K. Ploog, and G. H. Dohler, Phys. Rev. B 43, R12094 (1991). 
[13] G. K. Rasulova, Yu. A. Efimov and V. N. Murzin, J. Appl. Phys. 82, 3381 (1997). 

[14] K. J. Lug, H. T. Grahn, and K. H. Ploog, Phys. Rev. B 57, R6838 (1998). 
[15] M. Rogozia, S.W. Teitsworth, H.T. Grahn, and K.H. Ploog, Phys. Rev. B 64, 041308(R) 
(2001). 

[16] M. Rogozia, S.W. Teitsworth, H.T. Grahn, and K.H. Ploog, Phys. Rev. B 65, 205303 (2002). 

[17] M. Rogozia, H. T. Grahn, S. W. Teitsworth and K. H. Ploog, Physica B 314, 427 (2002). 

[18] A. Amann, A. Wacker, L. L. Bonilla, and E. Scholl, Phys. Rev. E, 63, 066207 (2001). 

[19] L. L. Bonilla, O. Sanchez and J. Soler, Phys. Rev. B 65, 195308 (2002). 

[20] L. L. Bonilla, in Nonlinear Dynamics and Pattern Formation in Semiconductors and Devices, 

edited by F.-J. Niedernostheide (Springer, Berlin, 1995), Chap. 1. 
[21] A. Wacker, in Theory of Transport Properties of Semiconductor Nanostructures, edited by E. 

Scholl (Chapman and Hall, London, 1998), Chap. 10. 
[22] A. Carpio, L. L. Bonilla, A. Wacker and E. SchoU, Phys. Rev. E, 61, 4866 (2000). 
[23] L. L. Bonilla and F. J. Higuera, SIAM J. Appl. Math. 55, 1625 (1995). L R. Cantalapiedra, L. 

L. Bonilla, M. J. Bergmann and S. W. Teitsworth, Phys. Rev. B 48, 12278 (1993); see Fig. 2. 

27 



